/* benchmark guru interfaces */
#ifndef OHMMS_FFT_BENCH_USING_FFTV0_H
#define OHMMS_FFT_BENCH_USING_FFTV0_H
#include <cmath>
#include <limits>
#include <fft/fft_v0.h>
#include <OhmmsPETE/OhmmsMatrix.h>
#include <Utilities/Timer.h>
using namespace std;

namespace APPNAMESPACE {

/** class to perform FFT benchmark
 * 
 * PREC precision of the data, double, float
 * T datatype: double, complex<double>, float, complex<float>
 * ffteng FFT engine as defined in FFT/fft.h
 * Using Matrix<T>
 */
template<typename PREC, typename T,bool domany, unsigned ffteng>
struct FFT_BENCH
{
  //typedef fft_v0::FFT_DFT<PREC,T,domany,ffteng> fft_engine_type;
  typedef FFT_DFT<PREC,T,domany,ffteng> fft_engine_type;
  ///save the enum
  enum {MANY_FFT=domany};
  ///data to be transformed
  Matrix<T> target;
  ///FFT engine which can do multiple ffts of domany==true
  fft_engine_type fft_engine;
  ///number of ffts
  int num_ffts;
  ///number of ffts
  int num_blocks;
  ///memory stride to perform blocked FFTs
  int stride_size;
  ///title of the benchmark
  string title;

  /** constructor
   * @param fft_dim FFT dimension
   * @param howmany number of ffts
   * @param buffersize fft_dim+2 for T=real
   * @param nb nunber of simultaneous FFTs 
   */
  inline FFT_BENCH(int fft_dim, int howmany, int buffersize, int nb=1)
    :target(howmany*nb,buffersize) 
     , fft_engine(fft_dim,target.data(),howmany)
     , num_ffts(howmany)
     , num_blocks(nb)
     , stride_size(buffersize*howmany)
  {
  }

  inline void setTitle(const string& a)
  {
    title=a;
  }

  inline void benchmark(int niters)
  {
    Timer myclock;

    myclock.restart();
    
    for(int iter=0; iter<niters; ++iter)
    {
      typename fft_engine_type::value_type * restrict ptr=target.data();
      for(int k=0; k<num_blocks; ++k, ptr+=stride_size)
        fft_engine.transform_forward(ptr,num_ffts);
    }
    fft_engine.timing[0]=myclock.elapsed()/static_cast<double>(niters);

    myclock.restart();
    for(int iter=0; iter<niters; ++iter)
    {
      typename fft_engine_type::value_type * restrict ptr=target.data();
      for(int k=0; k<num_blocks; ++k, ptr+=stride_size)
        fft_engine.transform_backward(ptr,num_ffts);
    }
    fft_engine.timing[1]=myclock.elapsed()/static_cast<double>(niters);
  }


  //inline void benchmark(int niters, double mflops)
  //{
  //  Timer myclock;

  //  myclock.restart();
  //  for(int iter=0; iter<niters; ++iter)
  //    fft_engine.transform_forward(target.data(),num_ffts);
  //  fft_engine.timing[0]=myclock.elapsed()/static_cast<double>(niters);

  //  myclock.restart();
  //  for(int iter=0; iter<niters; ++iter)
  //    fft_engine.transform_backward(target.data(),num_ffts);
  //  fft_engine.timing[1]=myclock.elapsed()/static_cast<double>(niters);

  //  ////fft_bench_many(fft_engine,target,niters,num_ffts);
  //  //report(mflops);
  //}

  inline double cpu_secs(int i) const { return fft_engine.timing[i];}

  void report(double mflops)
  {
    cout << fft_engine.fft_size << " " << num_ffts
      << " " << mflops/fft_engine.timing[0] << " " << fft_engine.timing[0] 
      << " " << title <<" for "  << endl;
      
    cout << fft_engine.fft_size << " " << num_ffts
      << " " << mflops/fft_engine.timing[1] << " " << fft_engine.timing[1] 
      << " " << title <<" back "  << endl;
  }

  inline void debug(bool print_details)
  {
    int fft_dim=fft_engine.fft_size;
    int howmany=fft_engine.num_ffts;

    double phase=2*M_PI/static_cast<double>(fft_dim);
    double norm=1.0/static_cast<double>(fft_dim);

    std::vector<double> sample(fft_dim);
    for(int i=0; i<sample.size(); ++i) 
      sample[i]=0.5*sin(phase*double(i))+0.1*sin(phase*2*double(i))+0.3*sin(phase*3*double(i))
        +0.5*sin(phase*4.1*double(i));;

    if(print_details)
    {
      cout.precision(12);
      std::cout << "  Input data size = " << fft_dim << endl;
      for(int i=0; i<fft_dim; ++i)
        std::cout << i << " " << sample[i] << endl;
    }

    for(int i=0; i<howmany; ++i)
      std::copy(sample.begin(),sample.end(),target[i]);

    fft_engine.transform_forward(target.data(),howmany);

    if(print_details)
    {
      std::cout << "  Forward-FFT data size = " << fft_dim << endl;
      const T* fftwed=target.data();
      for(int i=0; i<fft_dim; ++i)
        std::cout << i << " " << *fftwed++ << endl;
    }

    fft_engine.transform_backward(target.data(),howmany);
    target *= norm;

    T eps=10*numeric_limits<T>::epsilon();
    for(int i=0; i<howmany; ++i)
    {
      for(int k=0; k<fft_dim; ++k)
      {
        //cout << target(i,k) << " " << sample[k] << endl;
        if(abs(target(i,k)-sample[k])>eps)
          cerr << "WRONG " << target(i,k)-sample[k] <<  endl;
      }
    }

    cout.precision(6);
  }
};

/** declare convenient functions for fft_bench
 * @param engFFTENG fft engine, e.g., FFT_DFT<double,true,FFTW_ENG>
 * @param d data container of type MAT
 * @param niter number of interations
 * @param howmany howmany ffts
 */
template<typename FFTENG, typename MAT>
inline void fft_bench_many(FFTENG& eng, MAT& d, int niter, int howmany)
{
  Timer myclock;
  double dtf=0, dti=0;
  for(int iter=0; iter<niter; ++iter)
  {
    myclock.restart();
    eng.transform_forward(d.data(),howmany);
    dtf+=myclock.elapsed();

    myclock.restart();
    eng.transform_backward(d.data(),howmany);
    dti+=myclock.elapsed();
  }
  eng.timing[0]=dtf/static_cast<double>(niter);
  eng.timing[1]=dti/static_cast<double>(niter);
}

}
#endif
